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Mapping the neuronal circuits is essential to understand brain function. Recent technological advancements 
have made it possible to acquire the brain atlas at single cell resolution. Digital reconstruction of the neural 
circuits down to this level across the whole brain would significantly facilitate brain studies. However, 
automatic reconstruction of the dense neural connections from microscopic image stUl remains a challenge. 
Here we developed a spherical-coordinate based variational model to reconstruct the shape of the cell body 
i.e. soma, as one of the procedures for this purpose. When intuitively processing the volumetric images in 
the spherical coordinate system, the reconstruction of somas with variational model is no longer sensitive to 
the interference of the complicated neuronal morphology, and could automatically and robustly achieve 
accurate soma shape regardless of the dense spatial distribution, and diversity in cell size, and morphology. 
We believe this method would speed drawing the neural circuits and boost brain studies. 

Accurately drawing and quantitatively describing the map of neuronal circuits play an important role in 
understanding how brain works' ^\ Digital reconstruction of neuronal morphology from image stacks is 
essential and indispensable for drawing the neuronal circuits'* ''. Basically the connectivity between 
neurons can be inferred by digitizing the morphology of somas and dendrites, and then detecting the structural 
connections'". Recent progresses in molecular fluorescence labeling'^ and Golgi staining technique" and high- 
throughput imaging techniques'''"^^ allow imaging a mouse brain at micron or even submicron resolution, which 
in principle provide huge image datasets with adequate details to reveal the architecture ranging from one single 
neuron to a comprehensive brain-wide neural circuits. However, digitizing the morphology of neurons from these 
huge volume digital datasets is extreme difficult and appears as a bottle neck at neuronal level connectivity studies. 
One of the major obstacles is that in high resolution images, neuron and its fibers (called process) are too dense to 
perform an automatic segmentation for digital reconstruction. Shape reconstruction of somas, as an important 
part of the digital reconstruction, still remains a challenge. And this mainly results from the complexity of 
neuronal images, namely the huge variation in soma distribution density, soma size as well as the morphology 
of dendrites, all of which make it very difficult to describe the neurons using a single pattern. 

There are many methods developed to segment and reconstruct the shape of cells effectively from the micro- 
scopic images, such as watershed based methods^'"^', active contour based methods'"" '"'', gradient flow tracking 
algorithm^^, concavity detection methods'"''^', voronoi diagrams methods''" ''^, and the graph methods'"' ''^, all of 
which are able to accurately segment and reconstruct the shape of cells, each for specific purpose. However, none 
of these methods involve the reconstruction of soma's shape under the severe influence of dense and thick 
neuronal trucks. For this application, rayburst sampling method'*'' '"' has been proposed as a good tool for shape 
reconstruction of single neuron, and it works well in the segmentation of somas, dendrites, and even spines. 
However, they are not applicable to reconstruct the shape of clustered somas. The major obstacle here is that rays 
cast from the center point of a soma fail to end at the right soma boundary, instead they tend to wrongly end at the 
contacted adjacent soma boundary as designed to be stopped by the intensity threshold. 

In contrast to rayburst sampling method'", active contour model'" '"' might avoid this problem. By minimizing 
the designed variational model, active contour model can make the initial boundary elements move to the 
boundary where the biggest differences of signal intensity appear. Also, it is required that the boundary formed 
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Figure 1 | Comparisons of the robustness of VM_SCS and active contour model to the initial boundary, (a) Using active contour model to reconstruct 
the shape of the separately and closely positioned somas, the reconstructions (blue curves or blue dots) were largely determined by the initial boundaries 
(black circles), (b) With the same data and initial boundaries (black circles) in (a), almost identical curves (read curves) can be reconstructed using 
VM_SCS. The green curved (a) or the straight arrow (b) denoted as the trajectory that a boundary element moved, (c) The rate of accurate reconstructions 
of VM„SCS and active contour model were compared when the initial boundaries were set to the circles with the random centers and radius ranging from 
1.2 \im to 6 [im. (d) The intensity of gradient vector flow as a function of boundary elements that formed the boundary. The boundary was labeled by the 
blue curve in the embedded image, (e) Low isotropy of gradient vector flow makes the initial boundary (black circle) move to the red curve, the green 
curve and finally the blue dot. 



by the moving elements is smooth, preventing some elements from 
arriving on the boundary of other somas. However, active contour 
model still experiences difficulties in processing some complex 
images. The reason, in most cases, is that the low isotropy of gradient 
vector flow (GVF) in the images causes all boundary elements mov- 
ing towards a point. GVF here is regarded as the smoothing gradient 
of the boundary map pointing towards the boundary. Restricting 
moving directions of boundary elements may eliminate the effect 
of low isotropy of GVF. Inspired by this idea, here we have developed 
a special variational model, termed as VM_SCS (variation model 
based on sphere coordinate system), to reconstruct neuronal soma 
shape free from the influence of thick trucks in dense neural circuits. 
In this method, we convert the volumetric image data from the 
original Cartesian coordinate system to the sphere coordinate sys- 
tem. With the advantage of the sphere coordinate system, for each 
boundary element, its moving direction can be fixed by a preset polar 
and azimuth angle. This strategy decreases the dependence of our 
model on the initial boundary conditions, and thus provides a high 
robustness to the setting of the initial boundary. In addition, our 
model possesses the advantages of active contour model"'^ '"'. So our 
model can also accurately reconstruct the shapes of the clustered 
neuronal somas, and eliminate the interference of the thick dendritic 
trucks on soma shape reconstruction. With VM_SCS, the shape of 
neuronal somas with different size is also effectively reconstructed. 

Results 

The robustness of VM_SCS to the initial boundary. As active con- 
tour method*" is often frustrated by the initial boundary selection, we 
first explored the robustness of VM_SCS to the initial boundary in 
the shape reconstruction. For better demonstration, two 2-D images 
both having 40 X 40 pixels were chosen for our analysis (Fig. la or b). 
One contains two separate somas, and the other contains two 



contiguous somas. Note that, in the analysis of 2-D images, 
VM_SCS reconstructs the object's shape by solving the variational 
model in the polar system, and the corresponding parameters were 
set to the same as that of 3-D images. This was also suitable for the 
following analysis of 2-D images. 

In Fig. la, the results for the two images reconstructed using the 
active contour model showed strong dependence on the initial 
boundary (black circle). For the unsuitable initial boundary, its 
boundary elements converged to a fixed point (blue point, Fig. la), 
which usually results from the lacking of constraint on the direction 
that boundary element moves towards. This is illustrated by using the 
movement of a boundary element that started on the initial bound- 
ary, moved towards and reached the edge of soma, and continued to 
move along the edge and finally reached a fixed position, as shown 
with the green curved arrows (Fig. la). With the same images and the 
initial boundaries used in Fig. la, the shapes reconstructed using 
VM_SCS were almost identical (Fig. lb), which indicated that 
VM_SCS have more robustness to the selection of the initial bound- 
ary. In contrast to active contour model, VM_SCS makes all the 
boundary elements move along their own rays (green straight arrow. 
Fig. lb), in other words, the direction that each boundary element 
moves along is determined by polar angle and azimuth angle in 
sphere coordinate system, and thus VM_SCS can fix the directions 
of these elements. 

Furthermore, we quantified the robustness of VM_SCS to the 
initial boundary. In the calculations, all the initial boundaries were 
set to circles with uniformly and randomly positioned centers across 
a given square area. The square was located at the center of soma with 
6 |im side length. The size of the initial circle radius ranged from 
1.2 |im to 6 nm. This setting can avoid some multiple somas lying in 
the initial boundary, according to our results that the minimum value 
of the soma radius is about 3.5 |j,m. For each fixed radius value, 200 
simulated initial circles with different centers were generated and 
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Figure 2 | Shape reconstruction for the neuronal somas with dense spatial distribution, (a) Max-intensity projection of image stacks and shape 
reconstruction for the neuronal somas. Some reconstructions of the closely positioned somas are given (labeled in yellow arrows in (a)) and enlarged in 
(b), (c) and (d) respectively. Scale bar 50 |j.m in (a). 



divided equally into two parts. Then they were analyzed using 
VM_SCS and active contour model respectively. To give a quantitat- 
ive estimation, a 15% or less relative error in areas between automat- 
ically and manually reconstructed shapes is considered as an accurate 
reconstruction. This standard is also used in our following analysis. 
We calculated the rate of accurate reconstruction of VM_SCS and the 
active contour model with respect to a fixed radius of the initial 
circles with different centers. VM_SCS rendered an almost 
unchanged accurate reconstruction rate imder a wide setting of ini- 
tial boundary conditions, whereas the results using active contour 
model were largely dependent on the initial boundary (Fig. Ic). From 
these results, we conclude that VM_SCS has far more robustness 
than active contour model in terms of the initial boundary settings. 

Next we will explain why the reconstructions of active contour 
model were strongly dependent on the setting of initial boundary. 
Generally, in active contour model, gradient vector flow (GVF) and 
the smoothing force act simultaneously on the given initial boundary 
and make it evolve to the real boundary. GVF, the smoothing gra- 
dient of the boundary map, is supposed to point towards the real 
boundary in an ideal case, i.e. keeps high isotropy for our analysis. 
However, in Fig. Id, there is an abrupt change in the normalized 
intensity of GVF on the blue curve boundary. This intensity inhomo- 
geneity leaded to low isotropy of GVF in the neuronal soma, and 
damaged the shape reconstruction. It was illustrated by the example 
shown in Fig. le. Some GVF with the strong intensity can point 
towards the boundary, while others have the irregular directions 
(the black arrows in Fig. le). In this case, GVF with low isotropy 
made the unsuitable initial boundary (black circle) evolve to the red 
curves, to a part of the real boundary (green curve Fig. le) and finally 
to the given point (blue point in Fig. le). 

Shape reconstruction for the neuronal somas with different 
spatial distribution. To evaluate the capability of VM_SCS in 
reconstructing the shape of the neuronal somas with different 
spatial distribution, experimental image stacks with the size of 500 
X 500 X 20 volume pixels were analyzed. The neurons shown in 
Fig. 2a have different spatial distribution ranging from sparse to 
dense. It can be seen from the results that while VM_SCS is 
effective in reconstructing sparsely distributed neuronal somas 
(green arrows in Fig. 2a), it also yields high-precision shape 
reconstruction results for contiguous neuronal somas. As shown in 



Fig. 2a (yellow arrows) and the enlarged images of Fig. 2b-d, largely 
contiguous somas can be accurately reconstructed. Please note that 
we directly ignored the somas near the edge of image stacks, and thus 
some bright regions in Fig. 2a were not segmented. 

Furthermore, we also used simulation datasets to quantify 
VM_SCS for this purpose at different noise levels. Each dataset con- 
tained a pair of 5 \xm radius spheres with central distance ranging 
from 8 [im to 1 1 |im. The central distance here was denoted as d to 
show whether these spheres were closely positioned. The size of 
volume pixels was set to 0.5 X 0.5 X 0.5 |im', and each dataset 
had 60 X 60 X 60 volume pixels. The datasets were generated as 
follows. Inside the spheres, the volume pixel's value was set to be a 
Poisson random number whose expectation was the sum of the 
signal intensity Ig and the background Ij,. Outside the spheres, the 
value was Poisson random number with expectation If,. The signal- 
to-noise rate (SNR) is defined as the ratio of Ig to the standard 
deviation of lo+h- hi was a fixed value of 100 for all simulation 
datasets. Nine typical datasets were presented in Fig. 3a. Fig. 3b 
and 3c show the statistical information of the reconstructed sphere 
shape from the datasets. In shape reconstruction of single sphere, for 
each ray that a boundary element moved along, we computed the 
distance between the original position (center position of the sphere) 
and the position that the boundary element finally reached. Then the 
average and variance values of these distances were calculated 
(Fig. 3b and 3c). By comparing them to the ideal reconstruction, in 
which the average and variance should be 5 |im and 0 ^m^ respect- 
ively, we conclude that VM_SCS can reconstruct the shapes of the 
severely closely positioned somas at the suitable level of SNR. The 
reconstructions were further quantified in Fig. 3d by using the pre- 
vious standard that a 15% or less relative error in cubage between 
automatically reconstructed and real spheres. For severely closely 
positioned somas (d=8 |J,m), a correct rate of 80% can be achieved 
if signal-to-noise was more than 2. If signal-to-noise ratio dropped to 
1.5, the shapes of those closely positioned somas {d=9 \xm) can be 
effectively reconstructed with a correct rate of 96%. From the above 
results, we conclude that VM_SCS can effectively reconstruct shapes 
of separately and closely positioned somas at low level of signal-to- 
noise. 

Shape reconstruction of neuronal somas with different size and 
complicated morphology. We also used experimental datasets to 
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Figure 3 | Reconstruction performance of VM_SCS on the simulation somas with different spatial distribution and at different level of SNR. (a) x-y 

projections of nine pairs of simulation spheres with fixed radius of 5 [im presented. Here the alphabet of d represents the distance between the 
center positions of a pair of spheres, (b) The average lengths of rays in shape reconstruction of single sphere are shown and denoted by a point. The length 
of a ray means the distance between center position of the sphere and the reconstructed boundary position on this ray. For a frxed simulation parameters 
(d, SNR), 50 pairs spheres generated and their shapes reconstructed. Correspondingly, 100 points were shown for each group of parameters {d, SNR). 
(c) The variance of lengths of rays are shown when the distance d and SNR changes, (d) The accurate reconstruction rate as a function of the distance 
between the central positions of a pair of spheres for different SNRs ranging from 1 to 4. 



show that VM_SCS can reconstruct the shape of the neuronal somas 
with different size and complicated morphology. As shown in Fig. 4a, 
some neurons have long and thick trucks, resulting in the 
complicated morphology; the difference in the size of neuronal 
somas obviously exists. In the shape reconstruction of these 
neuronal somas, VM_SCS effectively eliminated the interference 
from those thick dendritic trucks and achieved accurate 
reconstruction (labeled as red surface in Fig. 4a). We also 
calculated the surface area and the cubage of digitized somas 
(Fig. 4b and c) respectively, using the protocol''\ The range of the 
surface area and the cubage of these somas are from 437 to 2373 |j,m^ 
and 680 to 8418 |j,m' respectively, which further verifies the large 
variance in the size of neuronal somas. 

Comparison of reconstruction results between VM_SCS and 
other methods. To further examine efficiency of VM_SCS, we also 
compared the reconstruction performance of active contour model, 
rayburst sampling and VM_SCS (Fig. 5). The experimental image 
has 1200 X 400 pixels with a pixel size of 1.2 X 1.2 |j,m^. In this 
image, many somas are closely positioned, and some of them are 
connected to the thick dendritic trucks, which is rather challenging 
for reconstruction. Before testing these methods, the positions of 
somas were manually given (black points in Fig. 5a-c), and the 
results were shown by red curves. In Fig. 5a, the shape of somas 
reconstructed using VM_SCS agreed well with those provided by 
the manual. By using the previous standard, among 129 somas 
there was only one soma whose shape can't be reconstructed. 
However, rayburst sampling method only obtained the boundary 
details of separately positioned somas, and cannot reconstruct the 
shapes of the neighboring somas contiguous to each other, as shown 
in the enlargement of the squares (Fig. 5b). 73 somas were 
successfully reconstructed (the correct rate of 57%) using rayburst 



sampling. When active contour model was used (Fig. 5c), due to the 
low isotropy of GVF in a neuronal soma caused by the disturbance of 
other somas or thick dendrite trucks, some reconstructions were 
points or open curves, as pointed out by the arrows in Fig. 5c. In 
this reconstruction, the correct rate of 77% was achieved. Obviously, 
we conclude that VM_SCS is more suitable for shape reconstruction 
for somas in complex neural images. 

Discussion 

In this paper, we have proposed the VM_SCS method to reconstruct 
the morphology of neuronal somas. In complicated image stacks of 
neurons, VM_SCS method can effectively reconstruct the shape of 
the clustered neuronal somas, and eliminate the interference from 
thick dendritic truck during reconstruction. We also show that 
VM_SCS is highly robust to the size variation of neuronal somas 
and the initial boundary in the reconstruction. 

Like the active contour model'*' '"*, our VM_SCS uses a combina- 
tion of the gradient signal and the smoothness force to drive the 
evolution of boundary elements to convergence. These convergent 
positions of boundary elements form the object's shape. It is clear, 
when some boundary elements evolve to the common boundary 
between two somas contiguous to each other or a soma and a thick 
dendritic truck, the gradient signal is weak, and the smoothness force 
is dominate, making these boundary elements stay at the common 
boundary, i.e., the correct positions. Based on this, our VM_SCS can 
reconstruct the shape of the clustered somas (Fig. 2, 3,) and eliminate 
the interference from the thick dendritic truck (Fig. 4). 

Different from active contour model''' *", our VM_SCS shows high 
robustness to the setting of the initial boundary. In active contour 
model'"*, the concept of gradient vector flow (GVF) is introduced to 
reconstruct the object's boundary. For applications with simple soma 
shape, the GVF diffuses in an isotropic manner, which keeps the 
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Figure 4 | Shape reconstruction for the somas of the neurons with different size and comphcated morphology, (a) Neurons with complicated dendritic 
trucks and the reconstructed soma shape (red). The density distribution of the surface area (b) and the cubage (c) of these reconstructed neuronal somas 
show the size diversity. 





Figure 5 | A comparison of three shape reconstruction methods on the experimental data. (a)The shapes reconstructed using VM_SCS under the 
condition of a given soma position, denoted by black points, (b) The shapes reconstructed using rayburst sampling that a rays were casted from a soma 
position (black points) and the threshold method were used to compute the length of the ray. (c) The shapes reconstructed using active contour model. All 
the initial boundaries were circles centered on their soma positions (black points) and whose radiuses were set to 5 |j.m for best results, (a) - (c) The results 
for two subareas (small black squares) selected for comparison, and further enlarged (big black squares). Note that, in the enlargement, the reconstructed 
shapes of different somas were distinguished using different colors, and the points or open curves that the arrows pointed towards still denote the 
reconstructed results. 
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Figure 6 | Scheme for the shape reconstruction routine of VM_SCS. (a) The sphere coordinate system with origin position O used to describe the signal 
from the image stacks, (b) The normalized signals in rays obtained via data sampling in sphere coordinate system. For each ray, it was equally divided into 
39 parts, and its signal is generated by averaging the signals of each dividing point (including endpoints) and its 6 neighbors, (c) The valid gradient 
signal of each ray was calculated from the curve in (b) with respect to the same ray, and generated by setting the negative value of the smoothing gradient 
signal to its converse value and setting the positive value to zero. The red curves in (b) and (c), corresponding to the red ray OP in (a), served as an example 
for illustrating how to obtain the gradient signal, (d) By iteratively computing the equation (7), which can minimize the energy function, the curve 
surface in sphere coordinate system was obtained from the signals in (c) . The third component value of a point in the curve surface is the distance between 
the computed boundary position to origin position O in the ray determined by the others components i.e., polar angle and azimuth angle, (e) Soma shape 
formed by the computed boundary positions in the rays. 



boundary elements move along their correct evolution directions. 
However, for applications with complicated neuronal morphology 
as described in this paper, this isotropic diffusion property of GVF 
might be broken for some initial boundaries, which leads to failed 
identification, as shown in the active contour model in Fig. 1 and 
Fig. 5c. While in VM_SCS method, the boundary element moves 
along its ray direction with the help of the sphere coordinate system, 
and naturally keeps the isotropic diffusion property during evolu- 
tion. This design contributes to high robustness of VM_SCS to the 
setting of initial boundary (Fig. 1 and Fig. 5a). 

Rayburst sampling'''' methods reconstruct the object shape 
according to the following steps: casting rays from the given point; 
computing the length the rays by a comparison of the signal intensity 
and a given threshold; using the length of rays for shape recovery. 
Due to a wide range of the number of sampling points in each ray, ray 
burst sampling can reconstruct the object's shape on different 
scales'*'"'"'. However, for clustered neuron in dense neural circuits, 
the signal intensity in the common boundary of the contiguous 
somas is stable, and does not experience a decrease below the thresh- 
old. So, rayburst sampling method cannot find the boundary of the 
contiguous somas (Fig. 5). VM_SCS obtains the optimal surface by 
minimizing the designed variational model (1st paragraph of 
Methods), shares the merits of active contour model (2nd paragraph 
of Discussion), and thus overcomes the inability of rayburst sampling 
in reconstructing the shape of clustered objects. 

The design of VM_SCS that the boundary element moves along its 
ray direction is based on the underlying fact that most neuronal 
somas have star-convex structure and can be presented in sphere 
coordinate system without information loss. This design contributes 
to the competitive results when using VM_SCS to analyze the com- 
plicated images (Fig. 1 and Fig. 5); meanwhile, this design also 
reduces the applied ranges of VM_SCS, namely, VM_SCS may be 
unable to completely reconstruct the non-star-convex shape. 
However, the case that a soma has non-star-convex shape is very 
rare according to our experience with thy-1 transgenic mouse; no 
exception is met across different brain areas including the cortex, the 
hippocampus, and epencephala. 

The robustness of VM_SCS against partial volume effects may be 
an important issue in shape reconstruction. However, in the shape 
reconstruction of cells field, due to the fact that molecular labeling 
technique characterized by specificity make partial volume effects 
negligible in imaging cells'^'" '", few algorithms have taken the issue 
of robustness against partial volume effects into account^'"'"'. 
Additionally, the issue involved in partial volume effects usually 
appears in the MRI and PET images^". Our VM_SCS aims at shape 
reconstruction of somas or other cells at present, and may not be 



directly applicable to MRI images because VM_SCS can't completely 
reconstruct some non-convex shapes. So, before applying VM_SCS 
to analyze MRI or PET images, the issue involved in partial volume 
effect should be considered in depth. 

For the high-throughput image stacks of neurons, the computing 
time of VM_SCS should be a matter of concern. In VM_SCS, sam- 
pling data in sphere coordinate system is necessary and wiU bring 
redundant data. However, compared to active contour model""*, 
VM_SCS fbces the evolution direction of boundary element and has 
a much faster convergence. Generally, using VM_SCS to reconstruct a 
neuronal soma takes about one second in an Inter(R)Xeon(R)CPU 
3.46 GHz computing platform. This speed should be similar to its in 
active contour model'"' according to our estimation. 

It should be noted that a rough estimation of position of neuronal 
soma is necessary before using VM_SCS for shape reconstruction. 
For neurons with complicated morphology, we have proposed 
NeuroGPS method^' to automatically locate the neurons across dif- 
ferent brain areas. We will combine VM_SCS with NeuroGPS 
method to automatically extract some features of neuronal soma 
including position, surface area, cubage, and so on. With the advan- 
tages of VM_SCS, we firmly believe that VM_SCS can also be applied 
to the reconstruction of morphology of complicated neuritis and 
spines, and should be a powerful computational tool for the digita- 
lization of neurons and drawing neural circuits. 

Methods 

Variational model in spherical coordinate system (VM_SCS) for shape 
reconstruction. In a spherical coordinate system with a fixed origin O, the position of 
a point P can be specified by the parameters r, (p, 6. Here, r, (p, 6 denote the distance of 
the point P and the point 0, polar angle, and azimuth angle respectively (See Fig. 6a). 
The morphology of neuronal soma can be modeled by the function r{(p,d) in the 
spherical coordinate system. The task of shape reconstruction of neuronal soma is 
transferred to estimating the function r{(p,9) from the image stacks. To characterize 
the function r{f,6), we assumed that the morphology of neuronal soma is smoothness 
to some extent, corresponding to smoothness of the function r{<p,6), i.e., a sharp 
change of r is impossible when tp, 6 vary slowly, and that the intensity of differential 
signal on the boundary of neuronal soma is stronger than that in the inside and 
outside of soma. Based on the above assumptions, by minimizing the following energy 
function 

E= 1^'" |^(=<i (1^) +=<2 (1^) )-/i£„,(r(p,e))|rfpd9 (1) 

we can obtain the function r{(p,6). Here the parameter a is the tradeoff between the 
smoothness of the function r{(p,0) and the loss of reconstructed details of neuronal 
soma. A big value of a will increase the smoothness of the function r{(p,0) and damage 
the reconstruction of the details of neuronal soma. Conversely, accurate 
reconstruction of the details is built on a sacrifice for the smoothness of the function 
r{(p,6)- The parameters czj, /5 are introduced to simplify the numerical calculation 
oir{(p,0). EextiKf'^)) is a function with respect to r{(p,Q), and can be written as 



SCIENTIFIC REPORTS | 4:4970 | DOI: 1 0. 1 038/srep04970 



6 



ds 



(2) 



where the function ^(5) is used to describe the intensity of the differential signal with 
respect to the position s in the ray determined by the parameters (p, 9. The detailed 
calculation of^(s) is in the section (see the entire processing of shape reconstruction, 
Step 2). Note that maximizing the function Eej^t{r{(p,9)) can be considered as finding 
the peak value of the function ^(5). Generally, the position 5 corresponding to the peak 
value of ^5) Ues on the boundary of neuronal soma. 

The numerical solution of VM_SCS. To minimize the energy function (1), we 
introduce the procedure that is similar to the active contour method*^. According to 
the well-known fact, if the curve surface r{(p,9) is the optimal solution that minimizes 
the energy function (1), it must satisfy the Euler equation 



2 V o(p^ c9 _ 



To solve the equation (3), we treat r{(p,9) as a function of time f, polar angle ^, and 
azimuth angle 6, and establish the new equation as 



E exp v(5) 

5El7([s„,,,-,^|) ^ ^ 

E exp(^^) 

set/(|s^,i.;]) ^ ^ 



Here, [] represents rounding the element to its nearest integer and f7(5} denotes the 
set including the position s and its 6-connected neighborhood, 1(5) is the signal value 
of position 5 from the image stacks. Each ray was uniformly divided into 40 sub 
intervals and its corresponding signal has been shown in Fig. 6b. The length of rays 
should be more than the biggest value of the radius of neuronal soma. 

Step 2. Calculation of the valid gradient signal in rays. For the given ray, the 
corresponding negative gradient signal was calculated from the signals v(5^) with m 
— 1, 2, M. The procedure is as follows 

a) Calculate the differential signal:Vv(s^) = v'(s„,-|-j) — v(5^„), m— 1,2,...,M. 

b) Using linear interpolation to calculate the value of Vv(5^), m — 1,1.5,2,. ..,M- 
1.5,M-1. 

c) Using the weighted averaging method to smooth the signalVv(5^) 



I VV(5^)| VV(S„) + 2Vv(s„_o.5) + 2Vv(5„ + o.5 ) 

|Vv(s„)|+4 



m = 1.5,2,2.5, •••,M-1.5 



dr((p,9,t) a. ( & 



r{f,9,t) d^r{ffi,t) 



(4) 



when r{(p,9,t) stabilizes, i.e., the right hand side of (4) approaches zeros, a solution of 
(3) can be achieved. By substituting equation (2) into equation (4), equation (4) can be 
rewritten as 



Vv(5^) is repeatedly computed until all small peaks inV\'(5f„) disappear. 
d) Calculate the negative gradient signal fromVv(5^) 

g(s„)= - min(Vv(s„),0), m = 1,2, • • • ,M- 1 



dr(wfi,t) ( d^r(m,e,t) ?^r{m,B,t) 

=<i ^3 +"-2 w 



dt 



3|^(s)exp(^- 



2(7^ 



(5) 



{s — r)ds 



To obtain the stable r{gj,6,t), we discretize the equation (5) 

At 



A(p/xi 



^nij,j ^i,j,k 



A9/0L2 



(6) 



where rjjj^ is an abbreviation of r{<pj, dj,t;^) and recognized as the function of the 
boundary element {(pj,9j), and 5„, is the sampling point in the ray determined by the 
parameters (pi, dj. Both A(p/a i and ^re set to one for simplifying equation (6), we 
can obtain 



"ij.k^i = nj.k - Ata(4r,j,t - rij^ i,t - rij_ i,t - r,_i - r,+i,j,k) 



2(72 



(Smj.j-n.i.k) 



(7) 



By iteratively computing until convergence, the numerical solution of r{(p,9) 

can be achieved. 

The entire processing of shape reconstruction. This section contains three steps. In 
step 1, we described signal sampling in the sphere coordinate system. From the 
sampled signals (Fig. 6b), we computed the valid nonnegative gradient signals 
(Fig. 6c) denoted as the function g in Eq. (7), which was illustrated in step 2. After 
obtaining the valid gradient signal, we iteratively computed Eq. (7) until convergence, 
and the surface reconstruction can be achieved (Fig. 6d and e), which was illustrated 
in step 3. Detailed description of these steps is as follows. 

Step 1. Signal sampling in the sphere coordinate system. In this section, the process 
of signal sampling in the sphere coordinate system with origin O is introduced. We 
discrete polar angle and azimuth angle respectively 

<Pi = '/Ni'' 0 = l,2,---,W,),e; = %2 7t a =1,2, ■■■Ml) 

with Ni and N2 being setting to 20 and 40 respectively. The position 5^_,j in the ray 
with polar angle (pi and azimuth angle 9j can be represented as 

Stnj.j — 0 + Vmi COS (pj COS 9j, COS ffij sin 9j, sin (p^j^ 

where with the sampling step of half of a voxel denotes the distance of the position 
Sff,,i,j-> and T denotes the transpose operation. The signal with respect to the position 
s^jj is computed by 



By these above operations, ^(5) is finally generated and regarded as the valid 
gradient signal (Fig. 6c). Considering that in our experimental image stacks, the size of 
volume pixels is 1.2 X 1.2 X 2.4 |j.m^ and sampling step is half of a volume pixels, we 
set M to 40 for all analysis, which is suitable for our experimental datasets. 

The smoothing operation in c) can achieve that the big and small value of the 
signalVv(5m,„) changes slowly and sharply respectively. This also means that we can 
keep the big peaks usually around the boundary of object shape, and discard the small 
peaks in the signal using a broad range of the number of iterative smoothness com- 
putation. In our dataset analysis, the iterative computation was carried out 100 times. 

In d), we select the negative gradient signal as the valid signal based on the fol- 
lowing facts: the signal intensity of intra-soma is more than that of extra-soma; the 
signal in sphere coordinate system was sampled from the intra-soma to the extra- 
soma. We can conclude that the differential signal around the boundary of neuronal 
soma should be negative, which indicates that the nonnegative differential signal are 
useless for the evolution of boundary element and its value can be set to zeroes. In 
addition, considering that the function (2) should be maximized for shape recon- 
struction, and that the position around the boundary of neuronal soma should 
approach the position of negative peak value of v(5^), the negative value of v{s^) was 
set to its opposite value. 

Step 3. Surface reconstruction. For the surface reconstruction, we repeatedly cal- 
culate the equation (7) untU all the boundary elements (^„ Oj) converge to the fixed 
positions r,j. These fixed positions formed the surface function r{<p,9) (see Fig. 6d and 
e). In equation (7), the parameter ^ was set to 



With (7 — 0.5, and Af is set to one. The setting of unknown parameters a and [}o was 
divided into two stages. In the first stage, to assure that the majority of boundary 
elements can reach the boundary of neuronal soma, the force of the smoothness in 
controlling evolution of boundary element should be reduced, and a and j^o were set 
to 0.2 and 0.8 respectively. In the second stage, for the minority of boundary elements 
that cannot achieve the boundary due to the small values of the nonnegative gradient 
signal, we enforce the smoothness for their evolution. The parameter a keeps 
unchanged and is set to 0. 

Brain tissue preparation. All experiments were performed in accordance with the 
guidelines of the Experimental Animal Ethics Committee at Huazhong University of 
Science and Technology. The THYl-EGFP-M line transgenic fluorescence mouse 
(adult (P41) male) was deeply anaesthetized, and trans car dially perfused with 50 mL 
phosphate buffered saline (PBS) (0.1 M) and 200 mL 4% paraformaldehyde (PEA). 
The whole brain was brought out, post-fixed in 4% PFA for 12 hours at 4 C, rinsed in 
PBS thrice (for 2 h, 2 h, and 12 h respectively), and then infiltrated in a graded series 
of GMA. Finally, the brain was embedded in gelatin capsules and polymerized for 
60 h at 60 'C. A detailed description of the samphng preparation can be found in Ref 
52. 

Imaging system. The whole brain was sliced and imaged using two-photon 
fluorescence micro-optical sectioning tomography system {2p-fM0ST)^^ that consists 
of dispersion compensation, fast scanning, microscope, and motion and sectioning 
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modules. The experimental lateral and axial spatial resolution of 2p-fMOST are about 

0. 45 ]im and 1.8 ^im respectively, and the size of volume pixel is down sampled to be 
1.2 X 1.2 X 2.4 [im^. These values ofthe parameters are sufficiently high to resolve the 
soma shape. The datasets were generated as follows. The image tiles were mosaicked 
to form a whole section. In each section, the inhomogeneous illustration patterns 
were removed using similar procedures proposed previously^^. An entire brain 
dataset contained about 6000 sections. A few parts selected from the entire dataset 
were used for our analysis. 
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